## ----setseedch8, echo = FALSE--------------------------------------------
set.seed(1982)

knitr::opts_chunk$set(
                      echo = FALSE,
                      results = "hide",
                      warning = FALSE,
                      message = FALSE,
                      error = FALSE
)


## ----loadlibsch8---------------------------------------------------------
library(reshape2)
library(dplyr)
library(tidyr)
library(lme4)
library(ggplot2)
library(scales)
library(boot)
library(rstan)
library(pander)
library(ggrepel)
library(stringr)
source("common_funcs.R")



## ----loaddatach8---------------------------------------------------------
if(file.exists("ch2_data.RData")){
  load("ch2_data.RData")
} else {
	warning("Not found")
}

## ----example-------------------------------------------------------------

## ----introexample--------------------------------------------------------
holder <- list()
for (i in 1:11) {
	holder[[i]] <- sheet1[,c("NEUTRAL", "CASEID","HANDDOWN",paste0("JUDGE",i),paste0("JUDGE",i,"DISS"))]	
	names(holder[[i]]) <- c("NEUTRAL", "CASEID","HANDDOWN","JUDGE","DISS")
}
holder <- do.call("rbind", holder)
holder <- subset(holder, !is.na(JUDGE))

## Decisions prior to the decision in MacDonald
holder <- subset(holder, HANDDOWN < as.Date("2011-07-06"))
nCases <- length(unique(holder$NEUTRAL))
holder$DISS <- ifelse(holder$DISS == "C", 1, 0)

holder <- holder %>%
    group_by(NEUTRAL) %>%
    mutate(w8 = 1 / length(unique(CASEID)))

### Create matrix
votemat <- acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                 value.var = "DISS")
votemat[is.na(votemat)] <- 0

holder$value <- 1
presencemat <- acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                     value.var = "value")
presencemat[is.na(presencemat)] <- 0

### Create weights
w8 <- apply(acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                  value.var = "w8"),
            1,
            function(x) unique(na.omit(x)))

### Sort the first two columns of a dataframe
### and remove duplicates
dedupe <- function(d) {
    pasteorder <- apply(d[,1:2], 1,
                        function(x) paste0(sort(x), collapse = "/"))
    ### Take the second one of the two
    d <- d[which(duplicated(pasteorder)),]
    d <- d[which(d[,1] != d[,2]),]
    return(d)
}
res <- list()
for (w in unique(holder$w8)) {
    v <- votemat[which(w8 == w),]
    p <- presencemat[which(w8 == w),]
    agree <- crossprod(v, v)
    sat <- crossprod(p, p)
    agree <- melt(agree,
                  varnames = c("judge1","judge2"),
                  value.name = "agree")
    sat <- melt(sat,
                varnames = c("judge1","judge2"),
                value.name = "sat")
    agree <- dedupe(agree)
    sat <- dedupe(sat)
    pos <- which(unique(holder$w8) == w)
    res[[pos]] <- merge(agree, sat, all = T)
    res[[pos]]$agree <- res[[pos]]$agree * w
    res[[pos]]$sat <- res[[pos]]$sat * w
}
res <- do.call("rbind", res)
example <- res %>%
    group_by(judge1, judge2) %>%
    summarize(agree = sum(agree),
              sat = sum(sat))
example <- subset(example, judge1 == 6 | judge2 == 6)
example$rate <- example$agree / example$sat
example <- example[order(example$rate),]



## ----getagreement--------------------------------------------------------
holder <- list()
for (i in 1:11) {
	holder[[i]] <- sheet1[,c("NEUTRAL", "CASEID","HANDDOWN",paste0("JUDGE",i),paste0("JUDGE",i,"DISS"))]	
	names(holder[[i]]) <- c("NEUTRAL", "CASEID","HANDDOWN","JUDGE","DISS")
}
holder <- do.call("rbind", holder)
holder <- subset(holder, !is.na(JUDGE))
holder$DISS <- ifelse(holder$DISS == "C", 1, 0)

holder <- holder %>%
    group_by(NEUTRAL) %>%
    mutate(w8 = 1 / length(unique(CASEID)))

### Create matrix
votemat <- acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                 value.var = "DISS")
votemat[is.na(votemat)] <- 0

holder$value <- 1
presencemat <- acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                     value.var = "value")
presencemat[is.na(presencemat)] <- 0

### Create weights
w8 <- apply(acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                  value.var = "w8"),
            1,
            function(x) unique(na.omit(x)))

### Sort the first two columns of a dataframe
### and remove duplicates
dedupe <- function(d) {
    pasteorder <- apply(d[,1:2], 1,
                        function(x) paste0(sort(x), collapse = "/"))
    ### Take the second one of the two
    d <- d[which(duplicated(pasteorder)),]
    d <- d[which(d[,1] != d[,2]),]
    return(d)
}
res <- list()
for (w in unique(holder$w8)) {
    v <- votemat[which(w8 == w),]
    p <- presencemat[which(w8 == w),]
    agree <- crossprod(v, v)
    sat <- crossprod(p, p)
    agree <- melt(agree,
                  varnames = c("judge1","judge2"),
                  value.name = "agree")
    sat <- melt(sat,
                varnames = c("judge1","judge2"),
                value.name = "sat")
    agree <- dedupe(agree)
    sat <- dedupe(sat)
    pos <- which(unique(holder$w8) == w)
    res[[pos]] <- merge(agree, sat, all = T)
    res[[pos]]$agree <- res[[pos]]$agree * w
    res[[pos]]$sat <- res[[pos]]$sat * w
}
res <- do.call("rbind", res)
res <- res %>%
    group_by(judge1, judge2) %>%
    summarize(agree = sum(agree),
              sat = sum(sat))
res$prop <- res$agree / res$sat
nJudges <- length(unique(holder$JUDGE))
nPairs <- sum(res$sat > 0)


## ----funnelplot, fig = TRUE, fig.cap = "Funnel plot of rates of agreement. The dotted lines indicate the range within which we would expect 95 percent of observations to fall if judges agreed at a constant rate", fig.width = 7, fig.height = 4----
plot.df <- subset(res, sat > 0)
plot.df$prop <- plot.df$prop * 100

### Get some intervals
### This is from https://gist.github.com/jkeirstead/7527432
calculate_funnel_limits <- function(n, theta, alpha=0.95, method="spiegelhalter") {

    ## Convert alpha into p
    p <- (1 - alpha)/2
    
    if (method=="exact") {
        ## Calculate the exact binomial counts
        ## Because this is a discrete distribution
        ## these limits will appear very jumpy
        upper <- qbinom(1 - p, 1:n, theta)/1:n
        lower <- qbinom(p, 1:n, theta)/1:n
    } else if (method=="spiegelhalter") {

        ## Calculate the confidence limits
        ## using Spiegelhalter's interpolation
        u1 <- qbinom(1-p, 1:n, theta)
        p.u1 <- pbinom(u1, 1:n, theta)
        u2 <- qbinom(1-p, 1:n, theta)-1
        p.u2 <- pbinom(u2, 1:n, theta)
        a <- ((1-p) - p.u2)/(p.u1-p.u2)
        upper <- (u2 + a*(u1-u2))/1:n

        l1 <- qbinom(p, 1:n, theta) - 1
        l1 <- replace(l1, l1<0, 0)
        p.l1 <- pbinom(l1, 1:n, theta)
        l2 <- qbinom(p, 1:n, theta)
        p.l2 <- pbinom(l2, 1:n, theta)
        a <- (p - p.l1)/(p.l2-p.l1)
        lower <- (l1 + a*(l2-l1))/1:n        
    }

    return(data.frame(x=1:n, up=upper, lo=lower))
}

require(reshape2)
tmp <- calculate_funnel_limits(n = max(plot.df$sat),
                               theta = weighted.mean(plot.df$prop, plot.df$sat) / 100,
                               method="exact")
tmp <- melt(tmp, id="x")
tmp$value <- tmp$value * 100

plot.df$label <- paste(num2judge(plot.df$judge1),
                       num2judge(plot.df$judge2),
                       sep="/")

plot.df$isSig <- FALSE
for (i in 1:nrow(plot.df)) {
    if (plot.df$sat[i] > 0) { 
        pval <- binom.test(x = ceiling(plot.df$agree[i]),
               n = plot.df$sat[i],
               p = weighted.mean(plot.df$prop,
                                 plot.df$sat)/100)$p.value
        plot.df$isSig[i] <- (pval < 0.04) ## for plotting
    }
    
}

p1 <- ggplot(plot.df, aes(x = sat, y = prop)) +
    scale_x_continuous( "Cases heard together") +
    scale_y_continuous("Rate of agreement (%)") +
    geom_hline(yintercept = weighted.mean(plot.df$prop, plot.df$sat)) + 
    geom_point(alpha = 0.3, position = position_jitter()) + 
    geom_line(data = filter(tmp, variable == "up"),
              aes(x = x, y = value),
              linetype = 2) +
    geom_line(data = filter(tmp, variable != "up"),
              aes(x = x, y = value),
              linetype = 2) +
    geom_text_repel(data = filter(plot.df, isSig &
                                  sat > 10),
                    aes(label = label),
                    size = 2.5) + 
    theme_uksc()

print(p1)


## ----exacttest-----------------------------------------------------------
res <- subset(res, sat > 0)
p <- weighted.mean(res$prop, res$sat)
x <- round(res$agree)
n <- round(res$sat)
good <- which(n > 0)
x <- x[good]
n <- n[good]
pvals <- sapply(1:length(n),
                function(i) {
                    binom.test(x = x[i], n = n[i], p = p)$p.value
                    })

alpha <- 0.05
sum(pvals < alpha)
mean(pvals < alpha)

bonferroni.adj <- alpha / sum(n > 0)
sum(pvals < bonferroni.adj)
mean(pvals < bonferroni.adj)
the_one <- which(pvals < bonferroni.adj)
the_one_cases <- res$sat[good][the_one]
the_one_prop <- round(100 * res$prop[good][the_one])


## ----irtplots, fig = TRUE, fig.cap = "Illustration of the parameters in item response theory models. ", fig.width = 7, fig.height = 12----
generic_irt_plot <- function(main = "", flip = FALSE, label = "") {
	plot(judgex, judgey,
		xlab = "",
		ylab = "Probability of a judge voting with the majority",
		axes = FALSE,
		type = "n",
		xlim = xlim,
             ylim = ylim,
             main = label)
	axis(2, at = c(0, 0.5, 1),
		labels = c("0%", "50%", "100%"),
		las = 1)
	abline(h = 0.5, col = 'grey')
	points(judgex, judgey,
		pch = 21,
		col = "black",
		bg = "grey")

	text(judgex, judgey, 
		labels = LETTERS[1:nPlotJudges],
		pos = 1,
		cex = 0.8)

        if (flip) {
            text(judgex, judgey + my.vadj, 
                 labels = paste0("(", round(judgex, 1), ")"),
                 pos = 1, 
                 cex = 0.8)

        } else {
            
	text(judgex, judgey - my.vadj, 
		labels = paste0("(", round(judgex, 1), ")"),
		pos = 1, 
             cex = 0.8)
        }
	### Mark out the third judge
	points(judgex[3], judgey[3],
		pch = 21,
		col = "black",
		bg = "darkred")
}

plot_pred_prob <- function(cp, discrim) {
	### Do the plotted probability
	theta <- seq(from = judgex[1], to = judgex[5],
		length.out = 100)
	inx <- discrim * theta - cp
	outy <- pnorm(inx)
	lines(theta, outy, lty = 2)

	### Add in a line to the cutpoint
	### 
	my.vadj2 <- diff(ylim) * 0.25
	### Draw the line

	arrows(x0 = cp/discrim, x1 = cp/discrim, 
		y0 = 0.5 - my.vadj2, y1 = 0.5,
		length = 0.1,
		lty = 1)
	text(x = cp/discrim, y = 0.5 - my.vadj2 - my.vadj,
		labels = bquote("Cutpoint: "~frac(.(cp),.(discrim))),
		cex = 1)

	### Add in lines showing the probability of C voting with majority
	## segments(x0 = judgex[3], x1 = judgex[3],
	## 	y0 = 0.5, y1 = pnorm(discrim * judgex[3] - cp))

	### And now from this point to a position to the right

	segments(x0 = judgex[3], y0 = pnorm(discrim * judgex[3] - cp),
                 x1 = judgex[3] + my.hadj, y1 = pnorm(discrim * judgex[3] - cp))

### Add some text
    xpos <- ifelse(discrim < 0, 0, -1.8)
    text(x = xpos, y = 0.9, labels = paste0("Location: ", cp), adj = 0, cex = 1)
    text(x = xpos, y = 0.9 - my.vadj, labels = paste0("Discrimination: ", discrim), adj = 0, cex = 1)
}

par(mfrow = c(4, 1), mar = c(1,3,1,2))
xlim <- c(-2, 2)
ylim <- c(-0.1, 1.1)
nPlotJudges <- 5
judgex <- seq(from = xlim[1] + diff(xlim)/8, 
	to = xlim[2] - diff(xlim)/8, 
	length.out = nPlotJudges)
judgey <- rep(0.5, nPlotJudges)
my.hadj <- 0.25
my.vadj <- diff(ylim) * 0.1
wrapn <- 50

### (a) Five judges, ordinary cut-point
### Choose a cut-point mid-way between 2 and 3
cp <- sum(judgex[2] + judgex[3])/2
### Choose the discrimination parameter
discrim <- 1

generic_irt_plot(main = paste0("Location: ", cp, "; Discrimination: ", discrim),
                 label = "(a)")
plot_pred_prob(cp = cp, discrim = discrim)

text(x = judgex[3] + my.hadj, y = pnorm(discrim * judgex[3] - cp),
	labels = str_wrap(paste0("When the cut-point is to the left of judge C, the probability of voting with the majority is high (", 
		round(100 * pnorm(discrim * judgex[3] - cp)),
		"%)"), wrapn),
	cex = 1,
	pos = 4)

### (b) Five judges, shift in cut-point
cp <- sum(judgex[1] + judgex[2])/2
discrim <- 1

generic_irt_plot(main = paste0("Location: ", cp, "; Discrimination: ", discrim), label = "(b)")
plot_pred_prob(cp = cp, discrim = discrim)

text(x = judgex[3] + my.hadj, y = pnorm(discrim * judgex[3] - cp) - my.vadj/2,
	labels = str_wrap(paste0("When the cut-point moves further to the left the probability of voting with the majority is higher (", 
		round(100 * pnorm(discrim * judgex[3] - cp)),
		"%)"), wrapn),
	cex = 1,
	pos = 4)

### (c) Five judges, positive discrimination parameter => more positive
discrim <- 2
cp <- sum(judgex[2] + judgex[3])

generic_irt_plot(main = paste0("Location: ", cp, "; Discrimination: ", discrim), label = "(c)")
plot_pred_prob(cp = cp, discrim = discrim)

text(x = judgex[3] + my.hadj, y = pnorm(discrim * judgex[3] - cp),
	labels = str_wrap(paste0("When the discrimination parameter becomes greater, distance to the cut-point has a greater effect on the probability of voting with the majority (", 
		round(100 * pnorm(discrim * judgex[3] - cp)),
		"%)"), wrapn), 
	cex = 1,
	pos = 4)

### (d) Five judges, negative discrimination parameter

cp <- sum(judgex[3] + judgex[4])/2
discrim <- -1

generic_irt_plot(main = paste0("Location: ", cp, "; Discrimination: ", discrim), flip = TRUE, label = "(d)")
plot_pred_prob(cp = cp, discrim = discrim)

text(x = judgex[3] + my.hadj, y = pnorm(discrim * judgex[3] - cp) - 0.09,
	labels = str_wrap(paste0("When the discrimination parameter becomes negative, judges further to the right are less likely to vote with the majority (", 
		round(100 * pnorm(discrim * judgex[3] - cp)),
		"%)"), wrapn),
	cex = 1,
	pos = 4)



## ----mappings, results = "asis"------------------------------------------
map <- read.csv("data/outcome-mapping.csv",
                check.names = FALSE)
pander(map,
       split.cell = 20,
       split.tables = Inf,
       caption = "Categorisation of left-wing outcomes according to area of law and identity of litigants (A) and (B)")


## ----irtprep-------------------------------------------------------------
holder <- list()
for (i in 1:11) {
	holder[[i]] <- sheet1[,c("NEUTRAL", "CASEID","HANDDOWN",paste0("JUDGE",i),paste0("JUDGE",i,"DISS"))]	
	names(holder[[i]]) <- c("NEUTRAL", "CASEID","HANDDOWN","JUDGE","DISS")
}
holder <- do.call("rbind", holder)
holder <- subset(holder, !is.na(JUDGE))
holder$DISS <- ifelse(holder$DISS == "C", 1, 0)

### Create matrix
votemat <- acast(holder, NEUTRAL + CASEID + HANDDOWN ~ JUDGE, 
                 value.var = "DISS")
### Which votes have zero variance when restricted
### to permanent judges (1-20)?
votevar <- apply(votemat[,1:20], 1, var, na.rm = TRUE)
votemat <- votemat[which(votevar > 0), 1:20]

stan.df <- melt(votemat,
                value.name = "y")
stan.df <- subset(stan.df, !is.na(y))
stan.df$caseid <- as.numeric(factor(stan.df$Var1))
stan.df$judgeid <- as.numeric(factor(stan.df$Var2))

N <- nrow(stan.df)
j <- stan.df$judgeid
k <- stan.df$caseid
J <- max(j)
K <- max(k)
y <- stan.df$y

## Create the X-Mat
doStatus <- function(app, resp, codes) {
     retval <- rep(0, length(app))
     case1 <- which(is.element(app, codes) & !is.element(resp, codes))
     case2 <- which(!is.element(app, codes) & is.element(resp, codes))
     retval[case1] <- 1
     retval[case2] <- -1
     return(retval)
}

sheet1$isCentralGovt <- doStatus(sheet1$APPTYPE,
                                 sheet1$RESPTYPE,
                                 c("E"))

sheet1$outcome2 <- case_when(grepl("A", sheet1$OUTCOME) ~ 1,
                             grepl("C", sheet1$OUTCOME) ~ -1,
                             TRUE ~ 0)

sheet1$pro_govt_ruling <- sheet1$isCentralGovt * sheet1$outcome2
###

sheet1$isPublicAuth <- doStatus(sheet1$APPTYPE,
                                 sheet1$RESPTYPE,
                                 c("D", "E", "F", "G", "H"))


sheet1$pro_publicauth_ruling <- sheet1$isPublicAuth * sheet1$outcome2
###

sheet1$id <- paste(sheet1$NEUTRAL,
                   sheet1$CASEID,
                   sheet1$HANDDOWN,
                   sep = "_")

aux <- data.frame(id = levels(factor(stan.df$Var1)))

aux <- merge(aux,
             sheet1,
             all.x = TRUE,
             all.y = FALSE)

left <- read.csv("left_outcomes.csv", stringsAsFactors = FALSE)
left <- unique(left[,c("CASEID", "LeftOutcomeReached")])

nrow(aux)
aux <- merge(aux, left,
             all.x = TRUE,
             all.y = FALSE)
nrow(aux)

hr <- read.csv("hr_cases_w_dissent.csv")
aux <- merge(aux, hr,
             by = "NEUTRAL",
             all.x = TRUE,
             all.y = FALSE)
aux$MajorityHR[is.na(aux$MajorityHR)] <- 0

aux <- aux %>%
    mutate(LeftOutcomeReached = case_when(LeftOutcomeReached == "Yes" ~ 1,
                                          LeftOutcomeReached == "No" ~ -1,
                                          TRUE ~ 0))

##
plaint <- read.csv("plaintiff_won.csv")

## Plaintiff wins = appwon * appellant = plaintiff
plaint <- plaint %>%
    mutate(AppIsPlaintiff = case_when(
               as.character(plaintiff) == as.character(APPNAME) ~ 1,
               as.character(plaintiff) == as.character(RESPNAME) ~ -1,
               TRUE ~ 0))

aux <- aux %>%
    mutate(AppWon = case_when(
               grepl("A", OUTCOME) ~ 1,
               grepl("C", OUTCOME) ~ -1,
               TRUE ~ 0))

aux <- merge(aux, plaint,
             by = c("NEUTRAL", "CASEID"),
             all.x = TRUE,
             all.y = FALSE)
aux$PlaintWon <- aux$AppIsPlaintiff * aux$AppWon
aux$PlaintWon[is.na(aux$PlaintWon)] <- 0

aux <- aux[match(levels(factor(stan.df$Var1)),
                 aux$id),]

stopifnot(aux$id == levels(factor(stan.df$Var1)))

mm <- model.matrix(~LeftOutcomeReached +
                       pro_govt_ruling +
                       pro_publicauth_ruling +
                       MajorityHR,
                   data = aux)

stan_data <- list(N = N, K = K,
                  J = J, j = j,
                  k = k, y = y,
                  X = mm, nPreds = ncol(mm))



## ----stanmodels----------------------------------------------------------
## This code taken from https://github.com/RobertMyles/IRT/tree/master/models/Stan
## and adapted

## What would be a decent prior
if (file.exists("cache/08_cached_results.RData")) {
    load("cache/08_cached_results.RData")
} else {
    f <- paste0(tempfile(), ".stan")
    stan_model <- 'data {
  int<lower=1> J; //Judges
  int<lower=1> K; //Cases
  int<lower=1> N; //no. of observations
  int<lower=1, upper=J> j[N]; //Legislator for observation n
  int<lower=1, upper=K> k[N]; //proposal for observation n
  int<lower=0, upper=1> y[N]; //vote of observation n
  int<lower=1> nPreds; // n. predictors
  matrix[K, nPreds] X; // model matrix
}
parameters {
  vector[K] alpha;
  vector[K] beta;
  vector[J] theta; // ideal point
  real<lower=0> sigma; // ideal pt variance
  real<lower=0> sigma2; // ideal pt variance
  vector[nPreds] kappa; // regression coefs
}
transformed parameters {
  vector[K] linpred; // discrimination
  linpred = X*kappa;
}

model {
  alpha ~ normal(0, 10);
  beta ~ normal(linpred, sigma2);
  kappa[1] ~ student_t(7, 0, 10);
  for (i in 2:nPreds) kappa[i] ~ student_t(7, 0, 2.5);
  theta ~ normal(0, sigma);
  theta[6] ~ normal(-1, 0.01); // constraint on Hale
  theta[7] ~ normal(1, 0.01); // constraint on Brown
  for (n in 1:N)
  y[n] ~ bernoulli_logit(theta[j[n]] * beta[k[n]] - alpha[k[n]]);
}

generated quantities {
  real<lower=0, upper=1> pcp;
  vector[N] pp;
  vector[N] cp;
  int<lower=0, upper=1> y_new[N];

  for (n in 1:N) {
      pp[n] = inv_logit(theta[j[n]] * beta[k[n]] - alpha[k[n]]);
      y_new[n] = bernoulli_rng(pp[n]);
      cp[n] = y_new[n] == y[n];
  }
  pcp = mean(cp[]);
}

'

    writeLines(stan_model, f)
    
    stan.fit <- stan(file = f,
                     data = stan_data, iter = 5000,
                     warmup = 2500, chains = 4,
                     pars = c("alpha", "beta",
                              "kappa", "theta",
                              "pcp", "sigma",
                              "y_new", "pp"),
                     thin = 5, init = "random",
                     verbose = TRUE, cores = 2,
                     seed = 1982)
   
    save(stan.fit, file = "cache/08_cached_results.RData")
}
s <- summary(stan.fit,
             probs = c(0.025, 1/12, 0.50, 11/12, 0.975))



## ----fitstatsch8, results = "asis"---------------------------------------
pcp <- s$summary[grep("pcp", rownames(s$summary)), "mean"]
pcp <- round(pcp * 100)
pcp.int <- s$summary[grep("pcp", rownames(s$summary)), ]
pcp.int <- pcp.int[c("2.5%", "97.5%")]
pcp.int <- round(pcp.int * 100)

null.pcp <- mean(stan_data$y)
null.pcp <- round(null.pcp * 100)

### GMP
getGMP <- function(predprobs, stan_data) {
    votes <- stan_data$y
    if (any(predprobs == 0)) {
        predprobs <- replace(predprobs,
                             which(predprobs == 0),
                             1e-16)
    }
    if (any(predprobs == 1)) {
        predprobs <- replace(predprobs,
                             which(predprobs == 1),
                             1 - 1e-16)
    }
    pijt_1<-log(predprobs)
    pijt_0<-log(1-predprobs)

    pijt_1<-pijt_1*(votes==1)
    pijt_0<-pijt_0*(votes==0)

    my.logLik<-sum(pijt_0,pijt_1,na.rm=T)

    my.A<-length(which(!is.na(votes)))

    GMP<-exp(my.logLik/my.A)
    GMP
}

my_predprobs <- extract(stan.fit, pars = "pp")$pp
gmps <- apply(my_predprobs, 1, function(x) getGMP(x, stan_data))
gmp <- round(mean(gmps), 2)
gmp.int <- round(quantile(gmps, c(0.025, 0.975)), 2)
## Alternately
null.gmp <- getGMP(rep(mean(stan_data$y), N), stan_data)
null.gmp <- round(null.gmp, 2)

### Produce a table
### PCP
### (Null PCP)
### Best predicted
### Worst predicted
### Proportion of cases that discriminate
### GMP

a <- data.frame(Model = "Item response model",
                GMP = gmp,
                `GMP 95% credible interval` = paste0("[",
                                                     paste0(gmp.int, collapse = ", "),
                                                     "]"),
                PCP = pcp,
                `PCP 95% credible interval`  = paste0("[",
                                                      paste(pcp.int, collapse = ", "),
                                                      "]"),
                check.names = FALSE)

b <- data.frame(Model = "Null model",
                GMP = null.gmp,
                `GMP 95% credible interval` = NA,
                PCP = null.pcp,
                `PCP 95% credible interval` = NA,
                check.names = FALSE)
fit.df <- rbind(a, b)

pander(t(fit.df), missing = "")



## ----idealplot, fig = TRUE, fig.cap = "Results of the ideal point analysis. Plotted points indicate best estimate; lines indicate 95 percent credible intervals. Positions of Hale and Brown constrained to be minus one and plus one respectively. See text for the source of the descriptions of judges as activist, moderate or restrained. ", fig.width = 6, fig.height = 10----

theta.bar <- s$summary[grep("theta", rownames(s$summary)), "mean"]
theta.lo <- s$summary[grep("theta", rownames(s$summary)), "2.5%"]
theta.llo <- s$summary[grep("theta", rownames(s$summary)), "8.333333%"]
theta.hi <- s$summary[grep("theta", rownames(s$summary)), "97.5%"]
theta.hhi <- s$summary[grep("theta", rownames(s$summary)), "91.66667%"]

theta.df <- aggregate(holder$DISS,
                         list(Judge = holder$JUDGE),
                         function(x) 1 - mean(x, na.rm = TRUE))
names(theta.df)[2] <- "dissent"
theta.df <- subset(theta.df, Judge %in% 1:20)

theta.df$theta <- theta.bar
theta.df$theta.lo <- theta.lo
theta.df$theta.hi <- theta.hi
theta.df$theta.llo <- theta.llo
theta.df$theta.hhi <- theta.hhi

theta.df$Judge.name <- num2judge(theta.df$Judge)
theta.df$Judge.name <- factor(theta.df$Judge.name,
                             ordered = TRUE)
theta.df$Judge.name <- reorder(theta.df$Judge.name,
                              theta.df$theta,
                              mean)

activist <- c("Hale", "Kerr")
restrained <- c("Sumption", "Clarke",
                "Hughes", "Toulson",
                "Reed", "Hodge")
moderate <- c("Neuberger", "Mance", "Wilson", "Carnwath")

theta.df$Type <- case_when(theta.df$Judge.name %in% activist ~ "Activist",
                           theta.df$Judge.name %in% restrained ~ "Restrained",
                           theta.df$Judge.name %in% moderate ~ "Moderate",
                           TRUE ~ "Unclassified")

p1 <- ggplot(theta.df, aes(x = Judge.name, y = theta,
                           shape = Type,
                           ymin = theta.lo, ymax = theta.hi)) +
    geom_linerange(size = inner.bar.size, alpha = 0.8) +
    geom_linerange(aes(ymin = theta.llo, ymax = theta.hhi),
                   size = outer.bar.size, alpha = 0.8) +
    geom_point(size = (inner.bar.size + outer.bar.size) * 2.5,
               stroke = 1.5) +
    scale_shape_manual("", values = my.pchs[c(1:2, 4, 3)]) +
    scale_x_discrete("") + 
    scale_y_continuous("Judge ideal points") + 
    coord_flip() +
    theme_uksc() +
    theme(legend.position = "bottom")


## Who is where?
print(p1)


## ----dicksoncomp---------------------------------------------------------
## Mean position of activist judges
theta <- extract(stan.fit, pars = "theta")$theta
activist.pos <- as.vector(unlist(theta[,num2judge(1:20) %in% activist]))
restrained.pos <- as.vector(unlist(theta[,num2judge(1:20) %in% restrained]))
moderate.pos <- as.vector(unlist(theta[,num2judge(1:20) %in% moderate]))

activist.mean <- round(mean(activist.pos), 2)
restrained.mean <- round(mean(restrained.pos), 2)
moderate.mean <- round(mean(moderate.pos), 2)

activist.int <- round(quantile(activist.pos, c(0.025, 0.975)), 2)
restrained.int <- round(quantile(restrained.pos, c(0.025, 0.975)), 2)
moderate.int <- round(quantile(moderate.pos, c(0.025, 0.975)), 2)



## ----arvindstirtoncomp---------------------------------------------------
arvst <- read.csv("arvind_stirton_estimates.csv", check.names = FALSE)
names(arvst)[2] <- "ArvindStirton"
theta.df <- merge(theta.df,
                  arvst,
                  by.x = "Judge.name",
                  by.y = "Judge.name",
                  all.x = TRUE,
                  all.y = FALSE)

arvst.cor <- round(cor(theta.df$ArvindStirton, theta.df$theta, use = "pairwise"), 2)



## ----cordissent----------------------------------------------------------
cordissent <- cor(theta.df$theta, theta.df$dissent,
                  use = "pairwise")
cordissent <- round(cordissent, 2)


## ----coefplot, fig = TRUE, fig.cap = "Coefficients of a model of case discrimination parameters. Thick lines indicate 83 percent confidence intervals. Thin lines indicate 95 percent confidence intervals. "----
kappas <- s$summary[grep("kappa", rownames(s$summary)), ]
kappas <- as.data.frame(kappas)
kappas$variable <- c("Intercept",
                     "Left outcome",
                     "Pro-central govt.",
                     "Pro-public authority",
                     "Pro-human rights")

kappas <- subset(kappas, variable != "Intercept")

kappas$isSig <- kappas[,"2.5%"] > 0 |
    kappas[,"97.5%"] < 0

kappas$bar.col <- ifelse(kappas$isSig,
                         "darkred",
                         "darkgrey")

ggplot(kappas, aes(x = variable, y = mean,
                   ymin = `2.5%`,
                   ymax = `97.5%`,
                   colour = bar.col)) +
    geom_linerange(size = inner.bar.size,
                   alpha = 0.8) +
    geom_linerange(aes(ymin = `8.333333%`,
                       ymax = `91.66667%`),
                       size = outer.bar.size,
                   alpha = 0.8) +
    geom_point(size = (inner.bar.size + outer.bar.size) * 2.5,
               shape = my.pchs[1],
               stroke = 1.5) +
    scale_x_discrete("Variable") +
    scale_color_identity() + 
    scale_y_continuous("Coefficient value") +
    geom_hline(aes(yintercept=0), lty=1) +
    coord_flip() +
    theme_uksc()


## ----illus, fig = TRUE, fig.width = 7, fig.height = 11, fig.cap = "Illustrative cases. Each row shows a case, judge locations, and 50 sampled probability curves. "----
## Which cases discriminate?
betas <- s$summary[grep("beta", rownames(s$summary)), ]
discriminates <- (betas[, "2.5%"] > 0 |
                  betas[, "97.5%"] < 0)

case.df <- data.frame(ID = levels(stan.df$Var1))
sheet1$nDiss <- apply(sheet1[,paste0("JUDGE",1:11, "DISS")],
                      1,
                      function(x) sum(x %in% c("A", "B")))
case.df <- case.df %>% 
    separate(ID, into = c("NEUTRAL", "CASEID", "HANDDOWN"),
             sep = "_") %>%
    mutate(HANDDOWN = as.Date(HANDDOWN)) %>% 
    left_join(sheet1 %>% dplyr::select(NEUTRAL, CASEID, HANDDOWN,
                                       Area, HR, nDiss))

case.df$discriminates <- discriminates
case.df$betabar <- betas[, "mean"]


illus <- case.df %>% arrange(desc(discriminates), desc(nDiss)) %>%
    distinct(NEUTRAL, Area, HR, nDiss, discriminates) %>%
    filter(nDiss > 1 & discriminates == TRUE)

contraindicators <- case.df %>% arrange(desc(discriminates), desc(nDiss)) %>%
    distinct(NEUTRAL, Area, HR, nDiss, discriminates) %>%
    filter(nDiss > 1 & discriminates == FALSE)

illus <- subset(holder, NEUTRAL %in% illus$NEUTRAL)
illus <- merge(illus, theta.df,
               by.x = "JUDGE",
               by.y = "Judge",
               all.x = TRUE,
               all.y = FALSE)

illus <- distinct(illus, CASEID, NEUTRAL, Judge.name, DISS, theta)

orig.illus <- illus
## Create an auxiliary data frame
## iteration by iteration, case by case
## the discrimination parameter and cutpoint
cp <- extract(stan.fit, pars = c("alpha", "beta"))
alpha <- melt(cp$alpha)
beta <- melt(cp$beta)
names(alpha) <- c("iter", "case", "alpha")
names(beta) <- c("iter", "case", "beta")
cp <- merge(alpha, beta, all = TRUE)
## Add on the neutral
cp$id <- levels(stan.df$Var1)[cp$case]
cp <- cp %>%
    separate(id, into = c("NEUTRAL", "CASEID", "HANDDOWN"),
             sep = "_") %>%
    group_by(NEUTRAL, iter) %>%
    summarize(alpha = mean(alpha), beta = mean(beta)) %>% 
    mutate(cp = alpha / beta) %>%
    filter(NEUTRAL %in% unique(illus$NEUTRAL)) %>%
    group_by(iter, NEUTRAL) %>%
    merge(data.frame(theta = seq(-2, 2, length.out = 100)),
              by = NULL)

cp <- cp %>%
    mutate(curve = pnorm(beta * theta - alpha)) %>%
    filter(curve > 0.05 & curve < 0.95) %>%
    distinct(NEUTRAL, iter, theta, curve) %>% 
    arrange(NEUTRAL, iter, theta)

illus$curve <- 0.5
tmp <- subset(illus, NEUTRAL %in% unique(illus$NEUTRAL)[1:3])
illus$NEUTRAL.pretty <- str_wrap(illus$NEUTRAL,
                                 7)
cp$NEUTRAL.pretty <- str_wrap(cp$NEUTRAL,
                              7)

## Create manual breaks
breaks_1 <- theta.df %>%
    arrange(theta) %>%
    filter(row_number() %% 2 == 0) %>%
    pull(theta)
labels_1 <- theta.df %>%
    arrange(theta) %>%
    filter(row_number() %% 2 == 0) %>%
    pull(Judge.name)
## Create a temporary data frame
breaks_2 <- theta.df %>%
    arrange(theta) %>%
    filter(row_number() %% 2 == 1) %>%
    pull(theta)
labels_2 <- theta.df %>%
    arrange(theta) %>%
    filter(row_number() %% 2 == 1) %>%
    pull(Judge.name)

illus_func <- function(df, aux) {
    judges_sitting <- unique(df$Judge.name)
    
    ggplot(df,
       aes(x = theta, y = curve)) +
    geom_line(data = subset(aux, iter < 51),
              aes(group = iter), alpha = 0.05) + 
    geom_point(aes(shape = factor(DISS),
                   colour = factor(DISS),
                   label = Judge.name),
               stroke = 1.5,
               size = (inner.bar.size + outer.bar.size) * 2.5) +
    facet_grid(NEUTRAL.pretty ~ .) + 
    scale_shape_manual(# guide = "none",
                       name = "",
                       values = c(4, 1),
                       labels = c("Dissented", "Concurred")) +
    scale_x_continuous("Judge ideal point",
                       breaks = breaks_1[labels_1 %in% judges_sitting],
                       labels = labels_1[labels_1 %in% judges_sitting],
                       sec.axis = sec_axis(~ .,
                                           breaks = breaks_2[labels_2 %in% judges_sitting],
                                           labels = labels_2[labels_2 %in% judges_sitting])) +
    scale_y_discrete("Case") + 
    ## geom_text(aes(label = Judge.name),
    ##           angle = 45,
    ##           hjust = 0,
    ##           nudge_y = 0.5) +
    scale_colour_manual(name = "",
                        values = c("darkred", "forestgreen"),
                        labels = c("Dissented", "Concurred"),
                        guide = "none") +
    coord_cartesian(xlim = c(-1.65, 1.25)) + 
    theme_uksc() +
    theme(legend.position = "bottom",
          strip.text.y = element_text(angle = 0),
          axis.text.x = element_text(angle = -90))

}

p5 <- illus_func(unique(illus), cp)



## ----hrillus, fig = TRUE, fig.cap = "Cases favouring a human rights interpretation", fig.height = 6, fig.width = 8----
hrillus_cases <- c("[2014] UKSC 38",
                   "[2010] UKSC 29",
                   "[2010] UKSC 33")

stopifnot(all(hrillus_cases %in% unique(illus$NEUTRAL)))

## Manually remove the one where 
hrillus <- illus_func(subset(illus, NEUTRAL %in% hrillus_cases),
                      subset(cp, NEUTRAL %in% hrillus_cases))
print(hrillus)

## Remove these cases from illus
illus <- subset(illus, !is.element(NEUTRAL, hrillus_cases))
cp <- subset(cp, !is.element(NEUTRAL, hrillus_cases))
             


## ----lrillus, fig = TRUE, fig.cap = "Cases favouring a left/right interpretation", fig.height = 6, fig.width = 8----
lrillus_cases <- c("[2011] UKSC 18",
                   "[2011] UKSC 58",
                   "[2014] UKSC 48",
                   "[2016] UKSC 38",
                   "[2011] UKSC 21",
                   "[2015] UKSC 2",
                   "[2012] UKSC 9",
                   "[2016] UKSC 35",
                   "[2016] UKSC 43")
stopifnot(all(lrillus_cases %in% unique(illus$NEUTRAL)))

lrillus <- illus_func(subset(illus, NEUTRAL %in% lrillus_cases),
                      subset(cp, NEUTRAL %in% lrillus_cases))
print(lrillus)

## Remove these cases from illus
illus <- subset(illus, !is.element(NEUTRAL, lrillus_cases))
cp <- subset(cp, !is.element(NEUTRAL, lrillus_cases))


## ----bothillus, fig = TRUE, fig.cap = "Cases favouring both interpretations", fig.height = 6, fig.width = 8----
bothillus_cases <- c("[2011] UKSC 23",
                     "[2015] UKSC 16",
                     "[2017] UKSC 41",
                     "[2015] UKSC 57")

stopifnot(all(bothillus_cases %in% unique(illus$NEUTRAL)))

bothillus <- illus_func(subset(illus, NEUTRAL %in% bothillus_cases),
                      subset(cp, NEUTRAL %in% bothillus_cases))

print(bothillus)

## Remove these cases from illus
illus <- subset(illus, !is.element(NEUTRAL, bothillus_cases))
cp <- subset(cp, !is.element(NEUTRAL, bothillus_cases))


## ----neitherillus, fig = TRUE, fig.cap = "Cases favouring neither interpretations", fig.height = 6, fig.width = 8----
neitherillus_cases <- c("[2009] UKSC 15 and [2009] UKSC 1",
                        "[2009] UKSC 15",
                        "[2010] UKSC 40",
                        "[2016] UKSC 4",
                        "[2017] UKSC 35")

stopifnot(all(neitherillus_cases %in% unique(illus$NEUTRAL)))

neitherillus <- illus_func(subset(illus, NEUTRAL %in% neitherillus_cases),
                      subset(cp, NEUTRAL %in% neitherillus_cases))

print(neitherillus)

## Remove these cases from illus
illus <- subset(illus, !is.element(NEUTRAL, neitherillus_cases))
cp <- subset(cp, !is.element(NEUTRAL, neitherillus_cases))

